function [M_sim, M_data,J_stat,Share_Sim,critical_value] = Jstat_end_myopic(betaStar,match)

NN = match.NN;
JJ = match.JJ;

N = match.N;
J = match.J;
M = match.M;
S = match.S;

% Parcel preferences
gamma = [1 betaStar(1,1:NN - 1)];
% Firm preferences
delta = [1 betaStar(1,NN:NN + JJ - 2)];
alpha = betaStar(1,NN + JJ - 1);

% Firm Value of parcels
ParcelUtility = match.ParcelIndex*gamma';
% Parcel value of firms (and leases)
FirmUtility = match.FirmMatrix*delta';


%%

tolHOLD = [];

SampleMMean_SimMoment       = zeros(M,JJ + NN,S);
SampleMVar_SimMoment        = zeros(M,JJ + NN,S);
SampleMCov_SimMoment        = zeros(M,JJ*NN,S);

SampleJMean_SimMoment       = zeros(J,NN,S);
SampleJCov_SimMoment        = zeros(J,NN*JJ,S);
SampleJVar_SimMoment        = zeros(J,NN,S);
SampleJDist_SimMoment       = zeros(J,NN*JJ,S);

SampleJntDist_SimMoment     = zeros(NN*JJ,S);
SampleCov_SimMoment         = zeros(NN*JJ,S);

Share_Sim                   = zeros(M*J,S);
ShareMean_Sim               = zeros(M,NN + JJ,S);
ShareFirm_Sim               = zeros(M,S);
ShareCov_Sim                = zeros(M,NN + JJ,S);

for s = 1:S

    shareIter = match.share0;
    shareKEEP = []; shareStKEEP = [];
    tol = 1; iter = 1;
    while tol > 1e-6 && iter < 300
        % Sub-market profits and index
        FirmUtility_sub = FirmUtility + alpha*shareIter + match.FirmError_Sim(:,s); 
        ParcelUtility_sub = ParcelUtility + match.ParcelError_Sim(:,s); 

        Aoffer = reshape(ParcelUtility_sub,[J,N/J]); 
        Baccept = reshape(FirmUtility_sub,[J,N/J]);

        % Sort preferences
        if match.out == 1
            rankA = SortingMatInC_out(Aoffer); 
            rankB = SortingResMatInC_out(Baccept);
        else
            rankA = SortingMatInC(Aoffer); 
            rankB = SortingResMatInC(Baccept);
        end

        % GS algorithm
        [~,FirmMatchOUT] = MatchAofferBaccept_GS_out(rankA,rankB,match.FirmMatchConstraint);
        [shareOUT,shareStOUT] = ShareMarket_st(FirmMatchOUT',match.MarketCount0,J);

        shareStKEEP = [shareStKEEP shareStOUT];
        shareKEEP   = [shareKEEP shareOUT];

        % tol = mean(abs(share0 - shareStOUT));
        if iter > 1
        %   tol = mean(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)))
            tol = sum(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)));
        else
        %   tol = mean(abs(share0 - shareStOUT))
            tol = sum(abs(shareIter - shareStOUT));
        end
        shareIter = shareStOUT;
%             shareIter = mean(shareStKEEP,2);

        iter = iter + 1;
    end
    tolHOLD = [tolHOLD tol];

    % Firm group royalty/clause means for moments
    [~,~,GpPMeanShare,GpFMeanShare,~,GpShareFirm,~,GpParcelCov,GpFirmCov] = GroupMoments_share2(match.ParcelIndex,match.FirmMatrix,FirmMatchOUT,J,match.MarketIndex0,match.MarketCount0,shareOUT);
    [~,~,PMean,FMean,JntDist,Cov,Pvar,Fvar,PFcovar,Pmean_firm,Pvar_firm,PFJntDist_firm,PFcovar_firm,~,~,~,~] = GroupMomentsSmall_new(match.ParcelIndex,match.FirmMatrix,FirmMatchOUT,J,match.MarketIndex0,match.MarketCount0,match.MomentSwitch);

    % Calculate the moments
    SampleMMean_SimMoment(:,:,s) = [FMean PMean];
    SampleMCov_SimMoment(:,:,s)  = PFcovar;
    SampleMVar_SimMoment(:,:,s)  = [Fvar Pvar];

    SampleJMean_SimMoment(:,:,s) = Pmean_firm;
    SampleJCov_SimMoment(:,:,s)  = PFcovar_firm;
    SampleJVar_SimMoment(:,:,s)  = Pvar_firm;
    SampleJDist_SimMoment(:,:,s) = PFJntDist_firm;

    SampleJntDist_SimMoment(:,s) = JntDist;
    SampleCov_SimMoment(:,s)     = Cov;

    % Calculate share - Market
    Share_Sim(:,s)          = shareOUT;
    ShareMean_Sim(:,:,s)    = [GpFMeanShare GpPMeanShare];
    ShareFirm_Sim(:,s)      = GpShareFirm;
    ShareCov_Sim(:,:,s)     = [GpFirmCov GpParcelCov];
end

%% 

JntDist_SimMoment    = sum(SampleJntDist_SimMoment,2)*(1/S);
Cov_SimMoment        = sum(SampleCov_SimMoment,2)*(1/S);

% Market Moments
MMean_SimMoment      = sum(SampleMMean_SimMoment,3)*(1/S);
MCov_SimMoment       = sum(SampleMCov_SimMoment,3)*(1/S);
MVar_SimMoment       = sum(SampleMVar_SimMoment,3)*(1/S);

% Firm Moments
JMean_SimMoment      = sum(SampleJMean_SimMoment,3)*(1/S);
JCov_SimMoment       = sum(SampleJCov_SimMoment,3)*(1/S);
JVar_SimMoment       = sum(SampleJVar_SimMoment,3)*(1/S);
JDist_SimMoment      = sum(SampleJDist_SimMoment,3)*(1/S);

% Market Structure Moments
Share_Sim           = sum(Share_Sim,2)*(1/S);
ShareMean_Sim       = sum(ShareMean_Sim,3)*(1/S);
ShareFirm_Sim       = sum(ShareFirm_Sim,2)*(1/S);
ShareCov_Sim        = sum(ShareCov_Sim,3)*(1/S);

M_sim = [JntDist_SimMoment; Cov_SimMoment; ...
    MMean_SimMoment(:); MCov_SimMoment(:); MVar_SimMoment(:); ...
    JMean_SimMoment(:); JCov_SimMoment(:); JVar_SimMoment(:); JDist_SimMoment(:)];
M = [match.JntDist0; match.Cov0; ...
    match.MMean_Moment; match.MCov_Moment; match.MVar_Moment; ...
    match.JMean_Moment; match.JCov_Moment; match.JVar_Moment; match.JDist_Moment];

% Full share set
M_share_sim = [Share_Sim; ShareMean_Sim(:); ShareFirm_Sim; ShareCov_Sim(:)];
M_share = [match.FirmShare0; match.GpMeanShare0; match.GpShareFirm0; match.GpCovShare0];

%%

disp 'J-stat'
J_stat = (N/J)*sum((1/(N/J))*([M_sim; M_share_sim] - [M; M_share]).*((1/(N/J))*([M_sim; M_share_sim] - [M; M_share])))

disp 'cannot reject the null that the moments fit if J is less than crit'
disp 'critical val'
critical_value = icdf('chi2',0.95,(size(M_sim,1) + size(M_share_sim,1)) - size(betaStar,2))

M_data = [M; M_share];
M_sim  = [M_sim; M_share_sim];

end